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Abstract. We present a formalism well adapted to the 
numerical study of the encounter of an ordinary main se- 
quence star with a massive black hole. Symmetry consid- 
erations, the use of a well adapted moving grid and a well 
adapted moving frame along with integration of the partial 
differential equations by means of pseudo-spectral meth- 
ods result in a very powerful and accurate tool. The hydro- 
dynamical equations are written in a moving frame which 
mimics the bulk of the movement of the fluid, resulting 
in very small relative velocities and a well suited spatial 
resolution throughout the calculations. Therefore, the nu- 
merical calculations are considerably simplified, smooth- 
ing in particular the Courant condition. Typical runs are 
performed within a few hours on a workstation with the 
high accuracy linked to the spectral methods. Predictions 
of the so-called affine are tested against this full numerical 
simulation. 
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1. Introduction 



The problem of the tidal influence of a massive black 
hole on an ordinary star is of great astrophysical inter- 
est, particularly in the case of close encounters. Carter 
and Luminet (1982, 1983) already suggested that the deep 
penetration of a star within the Roche radius of a black 
hole should strongly perturb its core and result in metal- 
enriched winds flowing out of the black hole, or even in 
helium detonation. Their model is a semi-analytical treat- 
ment based on the assumption of a linearized Lagrangian 
motion of the fluid. As a consequence, the shape of the 
star remains ellipsoidal. This model is referred to as the 
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affine model. It was numerically exploited in Luminet & 
Carter (1986) to predict the central density, temperature 
and entropy increases and in Luminet & Pichon (1989a, b) 
to estimate the additional nuclear energy release and the 
corresponding production of heavy elements. They showed 
in particular that, contrary to previous expectations, he- 
lium detonation by the triple-alpha proccess could prob- 
ably not occur although proton and alpha captures could 
change significantly the chemical composition of the star. 

A new motive for interest in this subject has been 
provided by recent theoretical considerations by Carter 
(1992) which show that the tidal disruption of an ordi- 
nary main sequence star is a conceivable scenario for the 
gamma ray bursts. Carter argues that the available en- 
ergy that might be radiated away through gamma rays, 
if a suitable transfert mechanism were available, would 
be of the order of the initial binding energy of the star, 
provided that the encounter is sufficently deep, with pene- 
tration factors of the order of 10. Such a mechanism might 
arise from an unstable shock formation due to deviations 
from affine behaviour. 

A qualitative description of the encounter is possible 
by making further approximations, depending on which 
part of the track the star is orbiting. One can split the 
movement of the star around the black hole into five dis- 
tinct phases. The first two phases correspond to fairly clear 
and reliable approximations. First, far from the black hole, 
the star is supposed to be in rough hydrostatic equilib- 
rium. When entering the Roche radius, the tidal acceler- 
ation tends to be dominant compared to self-gravity. The 
next three phases are the bounce, where pressure terms 
take over, then a rebound, which is an expansion with only 
tidal acceleration, and finally an ejection of material, pos- 
sibly driven by nuclear energy release. During these last 
three phases, the affine model becomes more and more 
unrealistic for several reasons: 

first, the geometry deviates strongly from an exact el- 
lipsoid when the star takes a double-wedged shape due to 
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the wringer effect of the tidal field (see e.g. Bicknell and 
Gingold 1983 or Carter 1992); 

second, hydrodynamical effects, like shocks, are likely 
to occur during the phase of very strong and rapid com- 
pression; 

third, if nuclear reactions are to be enhanced by the 
high temperatures achieved, or if a significant fraction 
of the total energy is released by electromagnetic waves 
(Carter 1992), the polytropic approximation becomes un- 
realistic. 

However, the affine model performs a very accurate de- 
scription of the movement of the star during the first two 
phases. 

To describe in more details the high compression 
phase, it then appears interesting to develop a hydrody- 
namical code which would be able to accurately follow the 
evolution of the star near the periastron and after, in or- 
der to give quantitative results on the possible generation 
of shocks and detonation of nuclear reactions. Several nu- 
merical investigations have tried to go beyong the affine 
treatment in working out the real deformation of the star 
near the black hole. The first investigation of this kind 
was made by Bicknell and Gingold (1983), using a 3-D 
Smoothed Particle Hydrodynamics (SPH) method. Their 
treatment was based on purely newtonian calculations. 
Their main result concerned the maximum heating and 
maximum compression of the star: they found less dra- 
matic effects than expected on the basis of the affine model 
and concluded that the triple-a could probably not det- 
onate although CNO reactions could change significantly 
the chemical composition. Morerecent SPH calculations 
were made by Evans & Kochanek, with a much better res- 
olution but in the axisymmetric approximation and not 
in full 3-D. Although their treatment is fully relativis- 
tic, they are interested only in the debris and not in the 
core itself, and therefore perform their calculation with 
a relatively small penetration factor. Further very recent 
publications address this problem again. Khokhlov et al. 
(1993a, b) report 3-D eulerian calculations where they ex- 
amine the energy and angular momentum tranfer to the 
star in order to check whether it might be disrupted or not. 
They do find a central density increase although quantita- 
tively different from the one predicted by Carter and Lu- 
minet (1983). However, they mention that their numerical 
method might not be reliable in every case due to its poor 
resolution. The same group (Frolov et al. 1994) followed 
up even more recently these calculations by including fully 
relativistic effects, both in the orbit and in the tidal field, 
using the analytic treatment by Marck (1983). They find 
quantitative differences from the non-relativistic case al- 
though the qualitative overall behaviour is similar. But 
the main effects described in those papers concern the 
outermost parts of the star, where the density is low, and 
the stripped mass does not exceed 10% of the total star 
mass. Laguna et al. (1993) have carried out 3-D relativis- 



tic SPH calculations to compute the evaporation of the 
star and the possible influence on emitted energy. They 
find significant deviations from the affine model which 
are compatible with the results of Bicknell and Gingold 
(1983) for the maximum compression. However, the SPH 
methods are known to be questionable (see e.g. Hernquist 
1993) since they effectively involve a high artificial viscos- 
ity that might lead to considerable entropy production in 
high compression phase and thus give erroneous results 
for the maximum central density which is crucial if one 
wants to determine what kind of nuclear reactions can be 
initiated or not. In a fully realistic description, one would 
expect energy dissipation in shocks, but it is still unclear 
how much of it will occur. It may be conjectured, at this 
stage, that the true outcome is intermediate between the 
predictions of the (strictly non dissipative) affine treat- 
ment and the (too highly dissipative) SPH treatment. 

The purpose of the present paper is to describe a new 
numerical approach which should help to settle these ques- 
tions. The essential idea is to combine the use of a mov- 
ing grid derived from the affine approach with the very 
powerful and reliable pseudo-spectral methods that have 
recently been developpcd for other purposes (Bonazzola & 
Marck 1990). The method we used is examplificd with one 
typical run. This is to be followed in the near future by an 
extensive study of the physical aspects of the encounter. 

A specific encouter may be characterized by a pene- 
tration factor (3. The definition of such a quantity is not 
unique. We choose the following form for [3 : 

P = R r /Rp (1) 

where 

R R = R*(M h /M*) 1/3 (2) 

is the Roche radius, Rp the periastron (minimum distance 
of the star to the black hole), R+ a characteristic radius 
of the star and Mh and M+ the mass of the hole and the 
mass of the star respectively. This definition agrees with 
the one taken by Laguna et al. (1993) but differs from 
that of Khokhlov et al. (1993a): their parameter ij can be 
identified with f3~ 3 / 2 . The details of the method are given 
in Sections 2 and 3. Preliminary results are presented in 
Section 4. Section 5 is devoted to discussion and prospects. 

2. Description of the method 

The study of tidal interactions needs a full 3-D calcula- 
tion. There are two mains categories of approaches: the 
Eulerian and the Lagrangian ones. In the first method, 
one writes the equations of motion with respect to a static 
frame. In the second method, one writes the equations of 
motion in a frame comoving with the matter. The Eule- 
rian method can be easily worked out but would require 
prohibitive number of mesh points to accurately describe 
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an inhomogeneous distribution of matter. The Lagrangian 
approach overcomes this drawback. However, in the case 
of multidimentionnal hydrodynamics, one has generally to 
take care of the formation of caustics. The Smoothed Par- 
ticles Hydrodynamics methods, which is a Lagrangian de- 
scription, can be applied whatever the matter distribution. 
However SPH is intrinsically highly dissipative and, hence, 
may lead to inaccurate results especially in the study of 
shock formation. 

We introduce an intermediate approach which com- 
bines the advantages of the Eulerian and Lagrangian 
methods. We solve the equations of motion in a mov- 
ing frame attached to the mean motion of the fluid. In 
the particular case of tidal interaction of a star and a 
massive black hole, the mean motion of the fluid is ac- 
curately described by the affine star model of Carter & 
Luminct (1983). After recalling the main features of the 
affine model, we write down the exact hydrodynamical 
equations in a general ellipsoidal coordinate system asso- 
ciated to its canonical frame. 

2.1. The affine model 

The fundamental hypothesis of the affine model described 
by Carter and Luminet (1983) is that the position of each 
cell of fluid with respect to a parallely propagated frame 
tied to the center of mass of the star can be described by 
a linear lagrangian transformation: 



(3) 



where r is the initial position vector of a fluid element and 
r is the current position vector in a frame which is paral- 
lely propagated along the orbit of the centre of mass of the 
star. Hence, the unknowns are the nine coefficients of 
the deformation matrix T>, which satisfy a system of sec- 
ond order differential equations which can be derived from 
a lagrangian. It turns out that, within this formalism, the 
star keeps ellipsoidal configurations. The principal axes 
of this ellipsoid are the eigenvectors of the matrix X>'2?. 
In the case of a planar orbit (i.e. newtonian approxima- 
tion or non-rotating black hole or an orbit lying in the 
equatorial plane of the Kerr black hole) the movement of 
the fluid in the z-direction and in the plane of the orbit 
decouple. As a consequence, T> has only 5 non-zero com- 
ponents. Solving the equations of motion for T> gives all 
the information about the physics of the encouter within 
the approximations made. The mass-density of each cell 
is given by 



P(r,t) 



P(r,0) 
\V\ 



and the velocity is 
v = T)r . 



(4) 



(5) 



In the particular case of a polytropic equation of state 
P ~ p 7 , the heat function of each cell obeys 



h{r,t) 



h(r,0) 
|D|7-i 



(6) 



When the star penetrates deeply inside the Roche radius, 
the surface of its equatorial section remains approxima- 
tively constant while the star undergoes a strong compres- 
sion in the z-direction. This induces an overall compres- 
sion. The maximum value reached by the central density 
is roughly given by: 



Po 



£2/(7-1) 



(7) 



The typical duration of this high compression phase is 



r m ~ /3"(7+i)/(7-i) 



(8) 



2.2. Hydrodynamics keeping the properties of the affine 
star model in mind 

We describe in this subsection the method we used to build 
our hydrodynamical code. Our point is to take advantage 
of most of the analytical results coming out of the affine 
model by writing the hydrodynamical equations in a well 
suited frame. We will be concerned in this paper with 
polytropic fluids, obeying P ~ p 7 . In that particular case, 
the enthalpy h is the right energy variable as shown in 
Marck and Bonazzola (1992) . The equations to be solved 
write: 

- mass conservation: 



d t p = -V- (pv) 

- energy conservation: 

d t h = -v-Vh-(rf- l)/iV • v 

- momentum conservation: 
d t v = -v-Vv — V(/i + $ + C) 

- Poisson equation 
A$ = AirGp 



(9) 



(10) 



(11) 



(12) 



where $ stands for the self-gravity potential and C for the 
tidal potential. Note that the continuity equation is not 
necessary when one uses an adiabatic relation P ~ p 7 . 

Let us now introduce the following coordinate trans- 
formation: 



X')-\0 Q-%t)J [x 



(13) 



whore V | y | are cartesian coordinates and where 
Q(t) is some 3x3 real regular matrix. Consider now the 
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Jacobian matrix J(t) associated with the previous trans- 
formation: 



(d T \ 

d x > 

8y> 

where 



1 J 



(d t \ 

d x 

\d z J 



(14) 



J 



1 
QQ- X X Q 



(15) 



and where a dot stands for time derivative. The time 
derivative operator is corrected for the grid movement, 
and thus is suitable for calculating relative velocity and 
acceleration with respect to the new frame. The relative 
velocity of the fluid with respect to the moving frame 
canonically associated to (r, X') reads 



v=Q- 1 v-Q- 1 QX / 



(16) 



and the hydrodynamical equations written in this new 
frame become 

- mass conservation: 



d rP =-V- (pv) - pd T log\Q\ 

- energy conservation: 

d T h = -v • Vh - (7 - l)h (V • v + d T log \Q 

- momentum conservation: 

d T v = -v • Vv - 2Q~ 1 Qv - Q^QX' 

- Q- 1 *G -1 V(/i + * + C) 



(17) 



(18) 



(19) 



where we have introduced V = 



By- 



It can be easily seen on the equations for p and h that 
the change of variables 



p(r'(t),t)=p(r'(t),t) \Q(t)\ 
and 

h(r'(t),t) = h(r>(t),t) IQ^r 1 



(20) 



(21) 



will absorb the extra terms in the equations [I?] and |l8|. 
The relevant variables we used are thus the components 
of the velocity relative to the moving grid v, the scaled 
density p and the scaled enthalpy h. The continuity and 
energy equations simplify further: 



d T p = -V- (pv) 



d T h = —v • Vh — (7 — 



(7 - l)h (V • v) 



(22) 



(23) 



and the Euler equation becomes: 

<9 r v = -v-Vv - 2Q~ 1 Qv - Qr x Q.X' 

Q- 1 'Q^V ' '" 



(24) 



$ + C 



We perform the calculations in a "pseudo-spherical" 
coordinate system linked to X' by the usual transforma- 
tions x' = r' smd' cos cf>' , y' — r' sin 8' sin </>' and z' = 
r' cos 9' because the relevant topology here is the one of a 
sphere. To be complete, we need to add the Poisson equa- 
tion, the equation of motion of the star and the explicit 
form of the tidal potential. Finally, to close the system 
of equations, we need to add a second order equation of 
motion for the 9 coefficients of the matrix Q. These last 
equations, whose choice is a priori arbitrary, will be set up 
in such a way that the grid motion is as close as possible 
to the mean matter motion. 

Many information can be drawn from the previous 
equations. The transformations we made were inspired of 
course by the affine model. Notice that if the motion of 
the fluid is exactly described by the equations of motion 
of the affine star model and if we give to the matrix Q 
the equation of motion of the deformation matrix T> (see 
section 2.1), the velocity of the matter with respect to an 
inertial frame linked to the centre of mass of the star is 
exactly given by 



Qr. 



(25) 



Hence, the vector r' is a constant and the relative speed 
v satisfies 

v = . (26) 

Moreover, p and h are then constant and the scalar fields 
p and h vary like powers of \Q\ : 

p(r'(t),t) = \Q(t)\- 1 p(r'(t = 0),0) (27) 
and 



h(r'(t),t) = \Q(t)\ L -^h(r'{t = 0),0) . 



(28) 



One possible choice for the equation of motion of Q would 
be to integrate simultaneoulsy the canonical affine equa- 
tions for the coefficients of the matrix T> and the hydro- 
dynamical equations. However, that way, the relative ve- 
locity of the fluid with respect to the moving grid at the 
surface of the star in the new frame would increase as the 
affine model becomes less and less realistic. We would not 
achieve a real improvement. We chose another way which 
has a number of great numerical advantages: we minimize 
the modulus of the velocity on the boundary of the star. 
The details of the procedure is given in the next Section. 
Doing so, we ensure that the star is approximately at rest 
during the first two phases of its track around the black 
hole, where its deformation is roughly ellipsoidal, and that 
the surface of the star is roughly given by r' = 1 in the 
moving frame. We thus save computation time and have 
better precision. 
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3. Numerical procedure 

Several tricks are used to reduce considerably the cost of 
the calculations. 

First, we take into account the symmetries of the prob- 
lem. The quadratic nature of all the potentials (tidal as 
well as gravitational) makes it possible to reduce the inte- 
gration domain. Taking advantage of the invariance with 
respect to the transformation: 

{x,y)->(-x,-y) (29) 

allows to make the calculation only on half a sphere. More- 
over, this problem is invariant under reflection with re- 
spect to the orbital plane 

z->-z. (30) 

The required integration domain is therefore only a fourth 
of a sphere. 

Second, the use of a well adapted moving grid has an 
important influence on the numerics. If the motion of the 
grid is chosen to be as close as possible to the mean motion 
of the fluid (which is always possible in the case of tidal 
interactions), the computed components of the velocity 
field almost vanish. As a consequence, the relative veloc- 
ity is very small, making the Courant condition (intrinsic 
to every hydrodynamical problem) very loose and hence 
allowing much larger time steps than otherwise required 
(by several orders of magnitude). 

Let us mention that, because the transformation giv- 
ing X'(t) from X is linear, the pseudo-singularities due 
to the use of the pseudo-spherical coordinate system 
(r' = 0,sin#' = 0) are of the same kind than the usual 
one. Therefore they can be handled like in Bonazzola and 
Marck (1990). 

As mentionned above, we use the pseudo-spectral 
method to solve our set of hydrodynamical equations. The 
method is a generalisation of the one described with ex- 
tensive details in Bonazzola and Marck (1990) which takes 
into account the symmetries (eqs. ^9| and [30]). The main 
advantage is that extremely high accuracy can be achieved 
without using any kind of artificial viscosity. Our integra- 
tion scheme is second order and semi-implicit in time. 

We used non dimensional variables and the non dimen- 
sioning procedure is as follows. Velocities are expressed 
in terms of the speed of light, parameters of the orbits 
in terms of the Roche radius, positions of fluid elements 
within the star in terms of the initial star radius. 

Let us give now a sketch of the integration algorithm. 
We first integrate explicitely all the terms of the right- 
hand side of eq. 25, except those containing . Then we 
calculate the qij in order to minimize the quantity: 

^boundary points (v ~ dtQ(t)" 1 Q(t)x') 

We then get a linear system, the solution of which 
gives the relevant . They are integrated in time with a 



second order scheme to get the q;j and the qy themselves. 
The initial conditions are provided by the canonical affine 
model. 

Finally we perform the implicit phase of the integra- 
tion. 

This gives the following scheme: 

v -> v J + dt (cl S v - c2 S^ 1 ) 

and 

h^h J +dt (cl S n -c2 S^ 1 ) 

where 

cl = (0.5dt J + dt J_1 ) /dt' 1 - 1 

c2 = 0.5dt J /dt J_1 

and S stands for the source terms without the grid accel- 
eration. 

Then, computation of the 'q\j by the method described 
before. Explicit integration of the remaining terms includ- 
ing q tj (only for v). 

v-> v -dt J Q _1 Qx' 

Finally, implicitation of the advection terms following the 
procedure: 

v J+i _ v 

— j = crd r v J+1 — crd r v 

at J 

and the same for v, where c is an adjustable parameter. 



4. A typical run 

4-.1. Run characteristics 

We present here the result of one specific run. This is to 
be followed in a next paper by an extensive study of the 
physics of the star core during the encounter, with partic- 
ular attention given to the possible development of shocks 
and detonation of nuclear reactions. It is not presently our 
aim to study the disruption and fate of the debris. 

Since the curvature around a super-massive black hole 
is small, and that the dimensions of a star are much 
smaller than the typical local curvature radius, a new- 
tonian treatment of the hydrodynamics is justified. How- 
ever, a relativistic treatment of the orbit and the tidal ten- 
sor is necessary for close encounters since they are known 
to have cumulative effects along the track (Luminet and 
Marck 1985). 

We take a polytropic star initially in spherical hydro- 
static equilibrium. We suppose that the orbit is parabolic 
and that the deformation of the star is exactly described 
by the model of Carter and Luminet down to the Roche 
radius. We thus integrate the canonical affine model from 
r=5 (where tidal effects are negligible and the spherical ap- 
proximation fairly good) to r=l. When the star crosses the 
Roche radius, we start the full hydrodynamical calcula- 
tion, with the initial conditions given by the affine model. 
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Fig. 1. Shape of the grid at different times. The grid is initially a fourth of a sphere. Each individual plot is centered on the 
corresponding location of the center of mass of the star on the track. Here = 1.5. One clearly sees the flattening in the z 
direction and the rotation of the principle axes in the orbital plane. 



The equation of state is polytropic in this first step, with 
7 = 5/3. The gravitational potential is not treated in the 
exact way by solving at each step the Poisson equation. 
We use instead a rough approximation which consists to 
keep the potential constant in time, its value being given 
by the initial model. We expect this approximation to be 
reasonable up to the point where autogravitation is defi- 
nitely negligeable, near the pericenter. 

The star is characterized by the density contrast be- 
tween the center and the boundary, which we chose to be 
10. The resolution of the calculation is 17 modes in the 
radial direction, 9 modes in 9 and 8 modes in <f>. For a 
specific encounter, we must specify a penetration factor 
0. The interesting range for is from a few to 10 or more, 
since, for such values, extremely high compressions are 



expected (Carter and Luminet 1983) and possible strong 
electromagnetic effects may occur (Carter 1992). 

The main difficulty in achieving numerical simulations 
of such encounters, is that the density contrast between 
the equilibrium configuration and the state of maximum 
compression may be as high as a hundred or even more, 
according to the prediction of the affine model: 

^ ~ 0^ = £3 if 7 = 5/3 (31) 
Po 

The code must resist this compression phase. 

The preliminary run we present here has = 1.5, cor- 
responding to rj <~ 0.54 in the definition of Khokhlov et 
al.. Taking such values for makes the approximation 
of a parabolic newtonian orbit very safe. But even this 
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relatively moderate penetration factor leads to a strong 
deformation of the grid, as illustrated on fig pi 



4-2. Comparison with the affine model 

The important parameter in the affine model is the central 
density. More precisely, the affine model makes a predic- 
tion on the maximum compression rate, (see eq. |3l| ) which 
can be compared in this hydrodynamical simulation with 
the ratio of the maximum to the initial central density. 
This is important if one wants to speculate on the pos- 
sible enhancement of nuclear reactions. Previous numeri- 
cal investigations of this problem by Bicknell and Gingold 
(1983) found a milder dependance on (3 than expected by 
the affine model, namely p m ax/po ~ ft 1 ' 5 , even for mod- 
erate values of (3. But at the same time, they find that 
a strong shock is formed and reverses the collapse of the 
stellar material toward the orbital plane. Of course, such 
a complex behaviour is not contained in the affine model 
and this might explain strong deviations for the value of 
the maximum compression. On the other hand, one could 
also argue that their method might not be reliable for 
very high compression and deformation, which can also 
lead to significant discrepancy. The numerical results by 
Khokhlov et al. ( 1993a, b) deal only with low values of 
(3 where essentially no compression occurs at all, except 
for one run, making comparisons difficult. Laguna et al. 
(1993), again using SPH, end up roughly with the same 
(3 dependance as Bicknell and Gingold (1983). Here we 
do not give any dependance of the compression rate with 
respect to the (3 parameter since we have only one run. A 
subsequent study including many runs and a careful study 
of the dependance on (3 of the results is still to come. How- 
ever, this first case seem to indicate a stronger compression 
than previously found in numerical calculations (see Fig. 

As expected, the use of a moving grid leads to much 
lower velocities than with a static grid, at least in the 
first stages of the encounter, as we could check by running 
the same calculation twice, once with a moving grid ac- 
cording to section 3 and once with a static grid. As the 
star approaches the periastron, it is very distorted and 
its internal motions do not coincide any more with what 
one can expect from an elliptical model. We present on 
figs ||, H and |^ the velocities relative to the moving grid, 
actually calculated by the code, at the periastron. Those 
velocities are the deviations to the pure affine behaviour. 
This presentation, although unusual, has the advantage 
of giving a feeling of the accuracy of the affine model. A 
final point we would like to make is that those relative 
velocities are plotted against the real radius, which gives 
the information about the geometry of star. The reader 
should bear in mind that the actual numerical calculation 
is made with an elliptical grid, so that in every direction 
the radius ranges from to 1. 




Fig. 2. Plot of the radial relative velocity (actually computed 
by the code) as a function of r for several values of the angles 
9 and <j> when the star reaches its periastron. 




Fig. 3. Plot of the relative velocity component vg (actually 
computed by the code) as a function of r for several values of 
the angles 9 and 4> when the star reaches its periastron. 

The time evolution of the matrix elements qij and their 
first and second derivatives as well as detQ are displayed 
on figs HI, @ and |. 

4-3. Comparison with other methods 

We already extensively discussed recent numerical results 
on the same problem. We would like to give in this sec- 
tion a few more technical elements, leading to quantitative 
comparison. 

The advantage of the SPH method is its versatility. 
Once the code is written, it can handle fairly easily almost 
any penetration factor. It seems also to be quick enough so 
that the number of particles is not a serious limitation. It 
seems possible within reasonable computing times to triple 
the number of particles. At least, the results of Laguna et 




Fig. 4. Plot of the relative velocity component (actually 
computed by the code) as a function of r for several values of 
the angles 9 and <f> when the star reaches its periastron. 



Matrix Q 




Fig. 5. Plot of the five components of the matrix Q with time. 
The overall behaviour is qualitatively identical to what one ex- 
pects from the affine model. Time is given in non-dimensional 
units, scaled by c/Rr 

al. do not seem to depend very much on the number of 
particles. Thus, one can have with this method a very 
good overview of the problem quickly. Moreover, for this 
particular problem, the study of the debris is well treated 
by SPH, even if qualitative results can be already given 
by a simple analysis of the geodesies. However, although 
it can in principle handle shocks as well, the precision 
is probably questionable, even if test cases are reproduced 
quite well and one has to be cautious with the quantitative 
results, especially in the high compression phase. 

The eulerian method used by Khokhlov et al. is more 
suited for distant encounters. It has the advantage of al- 
lowing a simple handling of boundary conditions and is 
in general less delicate than the spectral methods. But it 



Fig. 6. Same as fig |B|, but for the matrix Q. 

Matrix Q" 
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Fig. 7. Same as fig |B|, but for the matrix Q. 

Det|Qj 
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Fig. 8. Plot of det Q with respect to time. This picture sketches 
the behaviour of the central density very well, due to eq ^ij 
The maximum compression rate occurs at the maximum of det 

Q 
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is much more time consuming and fails to treat high-/? 
encouters accurately because the strong compression re- 
quires a prohibitive resolution. 

Our approach is much more accurate than any of the 
others. Boundary conditions and numerical implementa- 
tion in general are relatively complex but are now well 
mastered. We are able, with a moderate resolution, to 
compute a quite close encouter. Thus, the computing time 
is much less than for eulerian methods, owing also to the 
proper choice of the grid and frame we made. The time re- 
quired for the calculation presented is typically 2.5 hours 
on a Silicon Graphics Indigo workstation. Moreover, nu- 
merical difficulties grow very slowly with increasing /3, con- 
trary to others. 

5. Discussion and conclusion 

We have presented a new numerical approach to the tidal 
interaction of a star and a massive black hole, solving the 
hydrodynamical equtions in a moving grid with spectral 
methods, which proves very powerful and reliable. Our 
moving grid allows us to maintain a constant resolution 
within the star, whatever the deformation may be. By tak- 
ing advantage of most of the known analytical results to 
reduce the computation time and complexity, we achieve 
fully three dimensional calculations within a reasonable 
computing time, lower than the eulerian methods. More- 
over, this method has no artificial viscosity and thus virtu- 
ally no numerical dissipation, contrary to all other known 
methods (Laguna et al. 1993, Khokhlov et al. 1993a, b). 
We ran several calculations with moderate value of the 
penetration factor, but still greater than 1. 

We expect that the forthcoming extensive study is 
likely to give accurate additional results on the evolution 
of a star in a strong tidal field. In particular, the detailed 
balance of energy transfer could be examined as well as 
the oscillations subsequent to maximum compression in 
the case of low (3 encounters. For closer encounters, the 
possible occurrence of shocks and the dependance of the 
maximum density with [3 will be checked. 

However, the disruption of the star and the fate of the 
debris are unlikely to be accurately treated by this code. 
SPH for example seems to be much more adapted, but 
after all, the crescent shape of the debris found by Laguna 
et al. seems to be quite easy to predict by a simple geodesic 
calculation. 

Finally, this method might be applicable to a number 
of other physical applications including the oscillations of 
rotating stars. 

Acknowledgements. We would like to thank Brandon Carter 
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ject and for fruitful discussions in the course of this work. 
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A. Appendix 

In this appendix, we give the explicit form of the equations derived in Section 2 in pseudo-spherical coordinates. Recall that the 
continuity equation is useless in the case of a polytropic star. We do not write here the energy equation because it has the usual 
form. We drop the " and primes in this appendix to make things simpler to read. The symbol * stands for all the potentials in 
the Euler equation, that is, with the notations of Section 2: 

* = iqf^ + * +c (A1) 

where $ itself stands for the gravitational potential and C for the tidal potential. This latter potential has the usual expression 
in terms of the tidal tensor dj, in the newtonian approximation: 

C — — CtjXi^c j (A2) 
and 

Ci i = ""IT" ( 3x i x i - ^^i) • ( A3 ) 

The following equations are basically the components on a pseudo-spherical frame of equation 25 in the text. We applied 
the coordinate transformation followed by a projection onto the right vector basis. 



d T v r = - v r d r v r - (v g /r)d g v r - v^,/(rsm6)d c j > Vr + (vg + vl)/r 

+ 2v r (sin 2 9 A — cos 2 #933 /<?33 ) +2 cos 9 sin 9vg ( A + q 33 / q 33 ) + 2 sin 9 B 

+ (sin 2 6E - cos 2 6/qli) d r * + cos (9 sin (9 (l/g 2 3 + E) <9 9 * + smOFd^ + (sin 2 9H - cos 2 9q 33 /q 33 ) r (A4) 



d T Vg 


+ 
+ 


VrdrVg — (vgfr)dgvg - v^j(r sin 9)d^vg —v\ cot 9/r- (v r vg)/r 
2cos#sin#iv (q 33 /q 33 — A) + 2vg (cos 2 6 A - sin 2 0q 33 /q 33 ) +2cos9v <t> B 

cos 9 sin 9 (l/q 33 + E) 9 r * + (Scos 2 9 - sin 2 9/q 33 ) 9„* + cos6»F9 * + cos 9 sin 9r (933/933 + H) 


(A5) 


d T v i 


+ 
+ 


VrdrVg - (ve/rfdeVcj, - v 4 ,/(rsm9)d rl> v rl> - (v r v^)/r - v^vg cot 9/r 
2C (sin 9v r + cos 9vg) + 2v^D 
F(sin6»a r * + cos6'9 e *) - Gd^ + sin 9Ir 




(A6) 


where we adopted the following shorthand notations: 






A = 


(cos 4> 


sin0(giig2i -911921 - 912922 + 9121722) + cos 2 0(gi2g2i - <?ii<?22) + sin 2 0(qi 2 g2i 


- 911922)) /S 


(A7) 


B = 


(cos 4 


<sin0(gii<j22 - gn<?22 + gi2g2i - 912921) + cos 2 0(912922 - 1712922) + sin 2 0(giig 2 i 


- 911921)) /s 


(A8) 


C = 


(cos 4 


sin0(gng22 - gn<?22 + <ji2g2i - gi2<72i) + cos 2 0(gng 2 i - gii<?2i) + sin 2 0(gi 2 g22 


- gi2<?22)) /S 


(A9) 


D = 


(cos 4 


isin(/>(giig2i -911921 - gi2<?22 + qviqii) + cos 2 0(gi2g2i - 911922) + sin 2 0(gi 2 g2i 


- 911922)) /s 


(A10) 


E = 


(2 cos 


</)sin</)(giigi2 + qixqii) - cos 2 4(<]i2 + Q22) - sin 2 4{<lii + Q21)) /S 2 




(All) 


F = 


(cos 4 


sin(/>(g 2 2 - q\ x - q 2 21 + gf 2 ) + (cos 2 4 - sin 2 0)(gngi2 + g2ig22)) /S 2 




(A12) 


G = 


(2 cos 


4sin4(q 1 iq 1 2 + q2\q22) + cos 2 4{lii + Q21) + sin 2 4(<]i2 + 922)) /S 2 




(A13) 


H = 


(cos 4 


)sin(/>(g'iig2i - gng'21 - 912922 + 9121722) + cos 2 0(9121721 - 911922) + sin 2 0(g'i 2 g2i 


- 911922)) /<5 


(A14) 


I = 


(cos 4 


sin0(g'i 2 g2i - gi2<?2i + g'ng22 - gng'22) + cos 2 <Xg'iig2i - gng'21) + sin 2 <Xg'i2g22 - 


- 912922)) /5 


(A15) 



and 

S = 911922 - 912921 ■ (A16) 
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